PUBLPOL 218X Assignment 1

Code for Chapters 1.1 to 1.8

health_full = read_csv("https://chronicdata.cdc.gov/api/views/6vp6-wxuq/rows.csv?accessType=DOWNLOAD")
health_ca = filter(health_full, StateAbbr == "CA")

pge_q1 = read_csv("PGE_2019_Q1_ElectricUsageByZip.csv")

write_csv(health_ca, "health_ca.csv")
saveRDS(health_ca, "health_ca.rds")
save(health_ca, pge_q1, file = "working_datasets.rda")
load("working_datasets.rda")
save.image("progress1.rda")
load("progress1.rda")

year = 2019
quarters = 1:4
type = "Electric"
pge_19_elec = NULL
for(quarter in quarters) {
  filename = paste0("PGE_", year, "_Q", quarter, "_", type, "UsageByZip.csv")
  print(filename)
  temp = read_csv(filename)
  pge_19_elec = rbind(pge_19_elec, temp)
  saveRDS(pge_19_elec, "pge_19_elec.rds")
}

## [1] "PGE_2019_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q2_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q3_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q4_ElectricUsageByZip.csv"

pge_filter = filter(pge_19_elec, CUSTOMERCLASS %in% c("Elec- Residential", "Elec- Commercial"))
pge_select = select(pge_filter, -c(YEAR, COMBINED, AVERAGEKWH))
pge_group = group_by(pge_select, MONTH, CUSTOMERCLASS)
pge_summarize = summarize(pge_group, TOTALKWH = sum(TOTALKWH, na.rm = T))
pge_wide = pivot_wider(pge_summarize, names_from = CUSTOMERCLASS, values_from = TOTALKWH)
pge_tidy = pivot_longer(pge_wide, c("Elec- Commercial", "Elec- Residential"), names_to = "CUSTOMERCLASS", values_to = "TOTALKWH")

pge_summarize = summarize(pge_group, TOTALKWH = sum(TOTALKWH, na.rm = T), TOTALCUSTOMERS = sum(TOTALCUSTOMERS, na.rm = T))
pge_mutate = mutate(pge_summarize, AVERAGEKWH = TOTALKWH/TOTALCUSTOMERS)

pge_final = pge_19_elec %>%
            filter(CUSTOMERCLASS %in%
                     c("Elec- Commercial", "Elec- Residential")) %>%
            select(!c(YEAR, COMBINED, AVERAGEKWH)) %>%
            group_by(MONTH, CUSTOMERCLASS) %>%
            summarize(TOTALKWH = sum(TOTALKWH, 
                                     na.rm = T),
                      TOTALCUSTOMERS = sum(TOTALCUSTOMERS, 
                                           na.rm = T)) %>%
            mutate(AVERAGEKWH = TOTALKWH/TOTALCUSTOMERS)

pge_chart = pge_final %>%
            ggplot() +
            geom_bar(aes(x = MONTH %>% factor(),
                         y = TOTALKWH,
                         fill = CUSTOMERCLASS),
                     stat = "identity",
                     position = "stack") +
            labs(x = "Month",
                 y = "kWh",
                 title = "PG&E Territory Monthly Electricity Usage, 2019",
                 fill = "Electricity Type")
pge_chart

pge_chart %>% ggplotly() %>%
              layout(xaxis = list(fixedrange = T),
                     yaxis = list(fixedrange = T)) %>%
              config(displayModeBar = F)

pge_chart

Code for Chapter 1.9

ca_counties = counties("CA", cb = T, progress_bar = F)
st_crs(ca_counties)

## Coordinate Reference System:
##   User input: NAD83 
##   wkt:
## GEOGCRS["NAD83",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4269]]

projection = "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=ft +no_defs"
ca_counties_transformed = ca_counties %>% 
                          st_transform(4326) %>% 
                          st_transform(26910) %>% 
                          st_transform(projection) %>% 
                          st_transform(st_crs(ca_counties))

ggplot(ca_counties) + geom_sf()

leaflet() %>% addTiles() %>%
              addPolygons(data = ca_counties %>% 
                                st_transform(4326)) %>%
              addMarkers(data = ca_counties %>%
                                st_centroid() %>%
                                st_transform(4326))

bay_county_names <-
  c(
    "Alameda",
    "Contra Costa",
    "Marin",
    "Napa",
    "San Francisco",
    "San Mateo",
    "Santa Clara",
    "Solano",
    "Sonoma"
  )
bay_counties = ca_counties %>% filter(NAME %in% bay_county_names)
ggplot(bay_counties) + geom_sf()

ca_cities = places("CA", cb = T, progress_bar = FALSE)
bay_cities = ca_cities[bay_counties, ]
bay_cities_within = ca_cities %>%
                    st_centroid() %>%
                    .[bay_counties, ] %>%
                    st_set_geometry(NULL) %>%
                    left_join(ca_cities %>% select(GEOID)) %>%
                    st_as_sf()
leaflet() %>% addTiles() %>%
              addPolygons(data = bay_counties %>%
                                 st_transform(4326),
                          fill = F,
                          weight = 2) %>%
              addPolygons(data = bay_cities %>%
                                 filter(!GEOID %in% bay_cities_within$GEOID) %>%
                                 st_transform(4326),
                          color = "red") %>%
              addPolygons(data = bay_cities_within %>%
                                 st_transform(4326),
                          color = "green")